The critical role of ultra-low-energy vibrations in the relaxation dynamics of molecular qubits

Improving the performance of molecular qubits is a fundamental milestone towards unleashing the power of molecular magnetism in the second quantum revolution. Taming spin relaxation and decoherence due to vibrations is crucial to reach this milestone, but this is hindered by our lack of understanding on the nature of vibrations and their coupling to spins. Here we propose a synergistic approach to study a prototypical molecular qubit. It combines inelastic X-ray scattering to measure phonon dispersions along the main symmetry directions of the crystal and spin dynamics simulations based on DFT. We show that the canonical Debye picture of lattice dynamics breaks down and that intra-molecular vibrations with very-low energies of 1-2 meV are largely responsible for spin relaxation up to ambient temperature. We identify the origin of these modes, thus providing a rationale for improving spin coherence. The power and flexibility of our approach open new avenues for the investigation of magnetic molecules with the potential of removing roadblocks toward their use in quantum devices.

path covers all the symmetry direction explored in the IXS experiment. The x-axis labelling highlight the symmetry points of the crystal Brillouin zone (see Supplementary Fig. 1-c).

Supplementary Note 2 -Computational workflow
• DFT optimization is carried out with CP2K by optimizing both the lattice cell and the atomic positions. • Phonons are calculated by finite differentiating the forces with a two-point step of 0.01 Å • A and g tensors are calculated for 2000 randomly distorted structures using the ORCA package. The computational detail is specified below. • This is then used to train machine learning models in Python. The models with the best RMSE  are then used to calculate the first and second derivatives. The first and second derivatives for the first and second neighbours of Vanadium are computed both using machine learning and using DFT and the results are compared to ensure the quality of the coefficients calculated by machine learning (Supplementary Fig. 12). • The obtained derivatives and phonon modes are then passed to MolForge to compute the decay of the magnetization along or the plane . This is then fitted with exponential functions to obtain the relaxation times.

Cell Parameters
Optimized structure X-ray structure  Table 4. Optimized vs X-ray geometry comparison. Comparison between the cell parameters and bond lengths of the periodic DFT optimized [VO(TPP)] unit cell and the X-ray unit cell.
All the neural networks used in this work are made with the Tensorflow library and Keras API. The simple model has an input layer with 234 nodes, which translates to 3N Cartesian coordinates with N being the number of atoms in the molecule, and an output layer with 9 nodes, corresponding to the Cartesian tensor components of A or g. We tested the model with different numbers of hidden layers and the number of nodes in each hidden layer to obtain the optimized model for learning A and g tensors.
Each model is trained on a training set of 1600 configurations of randomly distorted [VO(TPP)] molecules. The random distortions are within the range of ±0.05 Å, applied to the equilibrium structure of [VO(TPP)]. The regularization hyperparameters are optimized for each model using a validation set of 158 configurations and then, each model's performance is evaluated with a test set of 200 configurations. DFT calculations are performed to obtain the needed data for the data set. A DKH-def2-TZVPP basis is used for V, N, and O atoms, and a DKH-def2-SVP basis is used for C and H atoms. All basis sets are decontracted and RIJCOSX is used as approximation for Coulomb and HF exchange. All DFT calculations are done using the hybrid functional PBE0 and a tight convergence in the ORCA software. Supplementary Figure 9 shows the learning curve of different models for the A and g tensors. Figure 9. Learning curves for the A and g tensors. The learning curve of machine learning models with different numbers of hidden layers for learning A tensor (a) and g tensor (b). The purple line and plus symbols report the learning curve of the model with 4 hidden layers; the green line and cross symbols report the learning curve of the model with 3 hidden layers; the blue line and asterisks report the learning curve of the model with 2 hidden layers; the orange line and unfilled squares report the learning curve of the model with 1 hidden layer; the yellow line and filled squares report the learning curve of the model with no hidden layer.

Supplementary
The models that performed the best for A and g tensors are chosen, with an error of 0.453 MHz for the A tensor and 8.236 × 10 −5 for the g tensor ( Supplementary Figs. 10-11).

Supplementary Note 3 -Spin relaxation simulations
The first-and second-order derivatives, and ′ , are calculated numerically from the A and g tensors computed with machine learning. A grid of 6×6 points is used to calculate each derivative and it is fitted using a 2D second-order polynomial expression with the linear least square method. The derivatives are then used to simulate the spin relaxation process of [VO(TPP)] and calculate the spin-phonon relaxation time 1 and 2 . We consider the quadratic term of the spin-phonon coupling Hamiltonian, which corresponds to a two-phonon process. The spin density ( ) of the system evolves in time as in which = ( − )/ℏ and the Redfield operator , is: The Dirac delta function ( ) is approximated using a Gaussian: The Redfield operator is converged to the harmonic limit by numerically taking the double limit The smearing in Eq. S4 and the q-point grid are varied to obtain the most accurate 1 and 2 profiles with respect to temperature.